using System;
using System.Collections;
using System.Collections.Generic;


namespace VerticalProcesses
{
    /// <summary>
    /// This is Georg's GRU which inherits from GRUClass...
    /// </summary>
    public class GRUGeorg: GRUClass
    {

        #region private variables
        
        // compartment members
        private Canopy RutterCanopy = new Canopy();
        private Soil HBVECSoil = new Soil();
        private SnowpackTwoLayer SnowPack = new SnowpackTwoLayer();
        
        // state variables

        private double _Mf;                     // potential snow melt rate, 
        private double _testV1;                 // test a variable for debugging
        private double _testV2;                 // test a variable for debugging
        private double _ccSurf;                 // cold content of snow surface layer, MJ m-2
        private double _ccPack;                 // cold content of snow pack layer, MJ m-2
        private double _tempSurf;               // temperature of snow surface layer, C
        private double _tempPack;               // temperature of snow pack layer, C


        // Model Parameters
        private double _Lat = 0.0;              // Latitude
        private double _Long = 0.0;             // Longitude
        private double _Timezone = 0.0;         // Timezone
        private double _CrownClosure = 1.0;     // Crown closure
        
        //snow related parameters
        private double _TT = 0.0;               // Temperature threshold for rain vs. snow
        private double _TTI = 1.0;              // Temperature interval for mixed rain and snow
        private double _MaxLiquid = 0.05;       // maximum liquid water fraction of snow, dimensionless
        private double _MaxSnowInt = 6;         // maximum snow interception
        private double _MaxRainInt = 4;         // maximum rain interception      
        private double _alphac = 0.0;           // canopy extinction coefficient
        private double _dia = 0.0;              // dia = tree diameter (m)
        private double _hs = 0.0;               // hs = stem height (m)
        private double _hc = 0.0;               // hc = height of the canopy (m)
        private double _n = 0.0;                // n= tree density (stems/(m^2))
        
        // soil related parameters
        private double _FC;                     //Field capacity of soil (mm)
        private double _BETA;                   //Exponent in relation between soil moisture and soil water release (dimensionless)
        private double _LP;                     //Soil moisture content as fraction of FC below which evaporation becomes supply-limited (dimensionless) 
        
        #endregion


        #region GRUGeorg Public Methods

        public override void Initialize(Hashtable InitialConditions)
        {
            throw new Exception("The method or operation is not implemented.");
            //initialize(InitialConditions);
        }

        public override void Update(DateTime TimeStep, double TimeStepHrs, Dictionary<string,double> MetForcing)
        {            
            
            double[] posClimate = (double[])METForcing["xyz"];
            
            //initialize flux variables
            waterRelease = 0.0;
            _Mf = 0.0;
            _testV1 = 0.0;
            _testV2 = 0.0;

            // climate forcings
            double gruTa = _Interpolation.LinearChangeWithElevation(posClimate[2],(double)METForcing["TA"], elevation[2], (double)METForcing["T_LAPSE"]);//-0.00546);//
            double gruP = _Interpolation.FractionalChangeWithElevation(1780.0, (double)METForcing["P"], elevation[2], (double)METForcing["P_LAPSE"]);//0.000317566);// ;                       
            double gruShortRad = ClearSkySolarRad(TimeStep, TimeStepHrs, (double)METForcing["T_LAPSE"]);           
            double ETpot = ETHamon(TimeStep, gruTa, gruShortRad, TimeStepHrs);
            double gruRain = gruP * rainFraction(gruTa);
            double gruSnow = gruP * (1 - rainFraction(gruTa));

            // update compartments
            // update canopy compartment
            UpdateCanopy(gruRain, gruSnow, _CrownClosure, ETpot,TimeStepHrs);
      
            // Calculate evaporation that reaches the soil or snow by subtracting the evaporation of  has already
            // been lost from the canopy from the potential evaporation  

            ETpot -= RutterCanopy.IntEvaporationRain - RutterCanopy.IntEvaporationSnow;
            
            // update snow compartment or direct throughfall directly to soil 
            if (RutterCanopy.ThrfSnow > 0.0 | SnowPack.Swe > 0.0)
            {

                if (SnowPack.Swe == 0.0)// new snow falling, zero CC and zero snowLiq assumed
                {
                    SnowPack.Swe = RutterCanopy.ThrfSnow;
                    SnowPack.TempSurf = 0.0;
                    SnowPack.TempPack = 0.0;
                    SnowPack.CcSurf = 0.0;
                    SnowPack.CcPack = 0.0;
                }

                //snowpack present
                // need to use this function even without snowfall to update cold content
                double snowPackEnergyInput = SnowPackEnergyInput(gruShortRad, TimeStepHrs, RutterCanopy.ThrfRain, gruTa);
                SnowPack.Update(RutterCanopy.ThrfRain, RutterCanopy.ThrfSnow, gruTa, snowPackEnergyInput);
                HBVECSoil.Infiltration = SnowPack.SnowOutflow;

            }

            else
            //no snow
            {
                HBVECSoil.Infiltration = RutterCanopy.ThrfRain;
            }

        
            // update soil compartment
            UpdateSoil(TimeStep, ETpot);
  
            //update water release from gru
            waterRelease =  HBVECSoil.WaterRelease;


            // save state variables
            
            //test variables here
            _testV1 = waterRelease;
            _testV2 = HBVECSoil.SoilMoisture;
            updateStateAndFlux();
        }


        public override void Change(Dictionary<string, double> Parameters)
        {
            parameters = Parameters;
            assignParameters();
            // add here more things that need to be done when parameters change
            // eg. empty canopy storage in case of a change in forest cover etc.

        }

        #endregion


        /// <summary>
        /// This is an example of a derived class from (GRUClass) which computes a
        /// one-dimensional water balance for a GRU. Note that we must override all 
        /// abstract methods.
        /// </summary>
        public GRUGeorg(string GRUID, double[] XYZ, double Area, double[] SlopeAspect, Dictionary<string, double> Parameters, bool output, string climateStation)
            : base (GRUID,XYZ,Area,SlopeAspect,Parameters, output, climateStation)
        {
            assignParameters();

            // initialize state variables (later add initial condidtions)

            updateStateAndFlux();
        }
       
        #region Private methods


        /// <summary>
        /// This function retrieves the parameters from the Parameters hashtable and assigns them to 
        /// 'real' variables
        /// </summary>
        private void assignParameters()
        {
            _Lat = (double)base.Parameters["latitude"];
            _Long = (double)base.Parameters["longitude"];
            _Timezone = (double)base.Parameters["timezone"];
            _CrownClosure = (double)base.Parameters["crownclosure"];
            _TT = (double)base.Parameters["tt"];
            _TTI = (double)base.Parameters["tti"];
            _MaxLiquid = (double)base.Parameters["whc"];
            _MaxSnowInt = (double)base.Parameters["maxsnowint"];
            _MaxRainInt = (double)base.Parameters["maxrainint"];
            _af = (double)base.Parameters["af"];
            _bf = (double)base.Parameters["bf"];
            _cf = (double)base.Parameters["cf"];
            _T0 = (double)base.Parameters["t0"];
            _alphac = (double)base.Parameters["alphac"];
            _dia = (double)base.Parameters["dia"];
            _hs = (double)base.Parameters["hs"];
            _hc = (double)base.Parameters["hc"];
            _n = (double)base.Parameters["n"];
            _FC = (double)base.Parameters["fc"];
            _BETA = (double)base.Parameters["beta"];
            _LP = (double)base.Parameters["lp"];
        }

        private void updateStateAndFlux()
        {
            fluxStateVars = new Hashtable();
            fluxStateVars.Add("swe", SnowPack.Swe);
            fluxStateVars.Add("coldcontent", SnowPack.CcSurf);
            fluxStateVars.Add("meltfactor", _Mf);
            fluxStateVars.Add("snowmelt", HBVECSoil.Infiltration);
            fluxStateVars.Add("intsnow", RutterCanopy.IntStorageSnow);
            fluxStateVars.Add("testV1", _testV1);
            fluxStateVars.Add("testV2", _testV2);
        }

        /// <summary>
        /// Computes the fraction of precipitation occurring as rainfall based on current
        /// air temperature 
        /// </summary>
        /// <param name="temperature">Curent air temperature</param>
        /// <returns>Fraction (0,1)</returns>
        private double rainFraction(double temperature)
        {
            if (temperature >= (_TT + _TTI)) return 1.0;
            if (temperature <= (_TT - _TTI)) return 0.0;

            // mixed rain/snow
            double frac = (temperature - _TT + _TTI) / (2 * _TTI);
            return frac;
        }


        private double ClearSkySolarRad(DateTime TimeStep, double TimeStepHrs, double tLapse)
        {
            // Time dependent Variables
            //  hour = hour of the day
            //  min = minute
            //  dn = day number  (1-365)
            //  lapse = temp lapse rate (K/m)

            // Site characteristics
            //  elev = elevation (m)
            //  Latitude = Latitude (degree)
            //  Latitude 	 = 49 + (20/60) + (40/3600)
            //  Longitude = longitude (degree)    For convention W is negative, E is positive
            //  Longitude	 = -115 + (44/60) + (56/3600)
            //  Lstandard = standard longitude (degree)
            //  Lstandard = -120
            //  aspect = aspect at site   (degree)
            //  slope = slope inclination at site (degree)

            // Constants
            // I0 = solar constant (W/m^2)
            // I0 = 1368
            // pa = mean atmospheric clear sky transmissivity
            // pa = .75
            // g = 9.80665 (gravity - m/s^2)
            // M = 0.0289644(molar mass of earth's air kg/mol)
            // R = 8.31432 (universal gas constant N*m/mol*K)
            // Tb = 288.15 (Standard temp K)

            // CALIBRATION COEFFICIENTS
            // alphac = extinction coefficient for canopy
            // Stand Characteristics
            //dia = tree diameter (m)
            //hs = stem height (m)  
            //hc = height of the canopy (m)
            //n= tree density (stems/(m^2))
            double Gf = 0.0;
            double hour = 0.0;

            for (int i = 0; i < TimeStepHrs; i++)
            {
                //average for every interval > 1 hour
                if (TimeStepHrs == 1) { hour = TimeStep.Hour; }
                else { hour = (double)i; }
                double min = 30.0;      //middle of time step
                double dn = TimeStep.DayOfYear;

                double cr = Math.PI / 180.0;                //convert inputs to radians
                double a = (180.0 - slope[1]) * cr;  //convert to south=0, east  = positive as required in Iqbal equations
                double beta = slope[0] * cr;
                double latRad = _Lat * cr;
                double timeZoneRad = _Timezone * cr;        //Long = Longitude*con
                double hourDec = hour + min / 60.0;         //change to decimal hour

                double gamma = 2.0 * Math.PI * dn / 365.25;   //day angle (rad)       (Iqbal 1983)            
                double dcl = (0.006918 - 0.399912 * Math.Cos(gamma) + 0.070257 * Math.Sin(gamma)    //declination (Iqbal 1983)
                    - 0.006758 * Math.Cos(2.0 * gamma) + 0.000907 * Math.Sin(2.0 * gamma)
                    - 0.002697 * Math.Cos(3 * gamma) + 0.00148 * Math.Sin(3 * gamma));
                double Et = (0.000075 + 0.001868 * Math.Cos(gamma) - 0.032077 * Math.Sin(gamma) //Et = equation of time (hour) (Iqbal 1983)
                    - 0.014615 * Math.Cos(2.0 * gamma) - 0.04089 * Math.Sin(2.0 * gamma)) * 229.18 / 60;
                double latime = hourDec + (_Timezone - _Long) / 15.0 + Et;    //local apparent time (hour) (Iqbal 1983) 15 results from 4/60....

                //hour angle (noon zero, morning positive)      (from John Richards model 2008)
                double omega = (12.0 - latime) * Math.PI / 12.0;

                //solar zenith angle        (Iqbal 1983)
                double costheta = Math.Sin(dcl) * Math.Sin(latRad) + Math.Cos(dcl) * Math.Cos(latRad) * Math.Cos(omega);
                if (costheta < 0.0) { Gf += 0.0; continue; } //return 0.0 when sun position is below horizon
                double theta = Math.Acos(costheta);

                //solar azimuth angle (east positive, south = 0, west negative)(Iqbal 1983), costheta = sin solar altitude
                double cospsi = (costheta * Math.Sin(latRad) - Math.Sin(dcl)) / (Math.Sin(theta) * Math.Cos(latRad));
                double psi = Math.Acos(cospsi);
                if (omega < 0.0) { psi = 0.0 - psi; } //east positive, south = 0, west negative(Iqbal 1983)

                //cosine of angle of incidence
                double ci = Math.Cos(beta) * Math.Cos(theta) + Math.Sin(beta) * Math.Sin(theta) * Math.Cos(psi - a); //Eq 1.6.5b Iqbal 1983 == Eq 9 Moore and Rowland 1992
                if (ci < 0.0) { Gf += 0.0; continue; } //avoid 0 division

                //As = area of shadow from a stem  (m^2)
                double As = _dia * _hs * Math.Sin(theta) / ci;     //Moore and Rowland 1992 Eq 8

                //Lc = path length through canopy (m)
                double Lc = _hc * Math.Cos(beta) / ci; //Moore and Rowland 1992 Eq 7


                //double tauCanopy = Math.Exp(-_alphac * _CrownClosure / Math.Cos(theta));

                //E0 = earth orbit eccentricity correction factor
                double E0 = 1.0 + 0.033 * Math.Cos(gamma - 0.048869);   //Iqbal 1983 and Grass Gis sun.r routine

                //barometric formula pp0=p/p0
                //double pp0 = Math.Pow((Constants.Tb / (Constants.Tb + tLapse * _XYZ[2])), ((Constants.g * Constants.M) / (Constants.R * tLapse)));
                //simpler formulation from r.sun in Grass:
                double pp0 = Math.Exp(-elevation[2] / 8434.5);
                //k0 = potential direct open site radiation (W/m^2) on slope
                double G0 = Constants.I0 * E0 * Math.Pow(Constants.pa, (pp0 / Math.Cos(theta))) * ci;  //Hock 1999

                //georg: Kf = potential k under forest (W/m^2), becomes kO at clearcuts since exp(0)=1
                Gf += G0 * Math.Exp(-_alphac * Lc) * Math.Exp(-_n * As); //Dan Moore personal communication              
            }
            return Math.Max(Gf / TimeStepHrs, 0.0);
        }
      
         
        /// <summary>
        /// Updates canopy storages and calculates interception
        /// and throughfall
        /// </summary>
        /// <param name="gruRain"></param>
        /// <param name="gruSnow"></param>
        /// <param name="_CrownClosure"></param>
        /// <param name="ETpot"></param>
        private void UpdateCanopy(double gruRain, double gruSnow, double _CrownClosure, double ETpot, double TimeStepHrs)
        {
           
                IntercRain(gruRain, _CrownClosure, ETpot);     
                IntercSnow(gruSnow, _CrownClosure, ETpot,TimeStepHrs);
        }
        
        /// <summary>
        /// Calculates snow interception and updates
        /// canopy compartment
        /// </summary>
        /// <param name="prec"></param>
        /// <param name="crownClosure"></param>
        /// <param name="etPot"></param>
        /// <returns>updates canopy</returns>
        private void IntercSnow(double prec, double crownClosure, double etPot, double TimeStepHrs)
        {
            double intRS;                                           // intercepted snow or snow
            double intEvap;                                         // interception evaporation
            etPot = 0.5 / 24.0 * TimeStepHrs;// 5.0 * etPot; //georg: a real fudge factor to increase sublimation...

            double catchPrec = prec * crownClosure * 0.7;           // snowcatch proportional to crownclosure change the 0.6 once real data!!!
            double maxInt = _MaxSnowInt * crownClosure;             // maximum interception storage for a GRU proportional to crown closure
            double intStTemp = RutterCanopy.IntStorageSnow + catchPrec - etPot;   //dry the canopy at ETpot

            if (intStTemp > 0.0) //canopy is wet throughout time step
            {
                intEvap = etPot; //let the potential evaporation evaporate
                if (intStTemp > maxInt) //canopy is at max interception storage
                {
                    intRS = intEvap + maxInt - RutterCanopy.IntStorageSnow;
                }
                else //still storage left in canopy
                {
                    intRS = catchPrec;
                }

            }
            else //canopy dries during time step or stays dry
            {
                intRS = catchPrec;
                intEvap = RutterCanopy.IntStorageSnow + catchPrec;
            }

            // update canopy snow properties
            RutterCanopy.IntStorageSnow = RutterCanopy.IntStorageSnow + intRS - intEvap;                  //canopy water balance
            RutterCanopy.ThrfSnow = prec - intRS;
            RutterCanopy.IntEvaporationSnow = intEvap;                                   // throughfall
        }

        
        /// <summary>
        /// Calculates rain interception and updates
        /// canopy compartment
        /// </summary>
        /// <param name="prec"></param>
        /// <param name="crownClosure"></param>
        /// <param name="etPot"></param>
        /// <returns>updates canopy</returns>
        private void IntercRain(double prec, double crownClosure, double etPot)
        {
            double intRS;                                           //intercepted rain or snow
            double intEvap;                                         // interception evaporation

            double catchPrec = prec * crownClosure * 0.5;           //raincatch proportional to crownclosure change the 0.6 once real data!!!
            double maxInt = _MaxRainInt * crownClosure;             //maximum interception storage for a GRU proportional to crown closure
            double intStTemp = RutterCanopy.IntStorageRain + catchPrec - etPot;   //dry the canopy at ETpot

            if (intStTemp > 0.0) //canopy is wet throughout time step
            {   
                intEvap = etPot; //let the potential evaporation evaporate
                if (intStTemp > maxInt) //canopy is at max interception storage
                {
                    intRS = intEvap + maxInt - RutterCanopy.IntStorageRain;
                }
                else //still storage left in canopy
                {  
                    intRS = catchPrec;
                }

            }
            else //canopy dries during time step or stays dry
            {       
                intRS = catchPrec;
                intEvap = RutterCanopy.IntStorageRain + catchPrec;
            }

            // update canopy rain properties
            RutterCanopy.IntStorageRain += intRS - intEvap;                  //canopy water balance
            RutterCanopy.ThrfRain = prec - intRS;
            RutterCanopy.IntEvaporationRain = intEvap;                                   // throughfall
        }
      

        /// <summary>
        /// Calculates potential evaporation according to Hamon (1961)
        /// retuns 0.0 during nighttime
        /// </summary>
        /// <param name="TimeStep"></param>
        /// <param name="gruTa"></param>
        /// <returns>potential evaporation</returns>
        private double ETHamon(DateTime TimeStep, double gruTa, double gruShortRad, double TimeStepHrs)
        {

            //if (gruShortRad == 0.0) // no evaporation during nighttime
            //    return 0.0;
            //// scale evaporation with temperature by calculating ETHamon 
            // as if for a day and prorate relative to daylength
            double P = Math.Asin(.39795 * Math.Cos(.2163108 + 2 * Math.Atan(.9671396 * Math.Tan(.00860 * (TimeStep.DayOfYear - 186)))));
            double dayLength = 24.0 - (24.0 / Math.PI) * Math.Acos((Math.Sin(0.8333 * Math.PI / 180.0) 
                + Math.Sin(_Lat * Math.PI / 180.0) * Math.Sin(P)) / (Math.Cos(_Lat * Math.PI / 180.0) * Math.Cos(P)));  //Approximate Day length
            //double SVD1 = 4.95 * Math.Exp(0.062 * gruTa);-->joels equation
            double SVD = (1324.0 / (gruTa + 273.2) * Math.Exp(17.27 * gruTa / (gruTa + 237.3)));//Dingman D-7, D8c
            return 0.55 * Math.Pow(dayLength / 12.0, 2) * SVD / 100.0 * 25.4 * TimeStepHrs / 24.0;
        }

        private void UpdateSoil(DateTime TimeStep, double ETpot)
        {
            double soilInfiltration = HBVECSoil.Infiltration;
            double evaporation = 0;
            double waterRelease = 0;
            double soilMoisture = HBVECSoil.SoilMoisture;

            if (soilInfiltration > 0.0)
            {

                double potEvaporationPerMM = ETpot / soilInfiltration;
                double infiltrationResidual = soilInfiltration - Math.Floor(soilInfiltration);
                double potEvaporationResidual = infiltrationResidual * potEvaporationPerMM;

                //Add 1mm (or residual) infiltration to soil at a time
                //Compute water release and actual evaporation
                double increment = 0.0;
                while (increment < soilInfiltration)
                {

                    //defaults assume increment < floor(soilInfiltration)
                    double actualInfiltration = 1.0;
                    double actualPotEvaporation = potEvaporationPerMM;

                    if (increment >= Math.Floor(soilInfiltration))
                    {
                        actualInfiltration = infiltrationResidual;
                        actualPotEvaporation = potEvaporationResidual;
                    }

                    double waterReleaseFraction = Math.Pow(soilMoisture / _FC, _BETA);

                    double thisWaterRelease = waterReleaseFraction * actualInfiltration;
                    actualInfiltration = (1.0 - waterReleaseFraction) * actualInfiltration;
                    //	waterRelease += waterReleaseFraction * actualInfiltration;
                    //  soilMoisture += (1.0-waterReleaseFraction) * actualInfiltration;

                    if (soilMoisture + actualInfiltration > _FC)
                    {
                        thisWaterRelease += soilMoisture + actualInfiltration - _FC;
                        actualInfiltration = _FC - soilMoisture;
                        //waterRelease += soilMoisture - FC;
                        //soilMoisture = FC;
                    }

                    waterRelease += thisWaterRelease;
                    soilMoisture += actualInfiltration;

                    //Evaporate
                    double thisEvaporation = CalculateSoilEvaporation(actualPotEvaporation, soilMoisture);
                    evaporation += thisEvaporation;
                    soilMoisture -= thisEvaporation;
                    increment += 1.0;

                }

            }
            else // no infiltration
            { 
                //evaporate
                double thisEvaporation = CalculateSoilEvaporation(ETpot,soilMoisture);
                evaporation += thisEvaporation;
                soilMoisture -= thisEvaporation;
            
            }
            // update actual evaporation and soil moisture arrays

            HBVECSoil.Evaporation = evaporation;
            HBVECSoil.SoilMoisture = soilMoisture;
            HBVECSoil.WaterRelease = waterRelease;
            
        }
        /// <summary>
        /// Calculates and returns soil evaporation
        /// </summary>
        /// <param name="potEvaporation"></param>
        /// <returns></returns>
        private double CalculateSoilEvaporation(double potEvaporation, double soilMoisture)
        {
            if (SnowPack.Swe > 0.0 && _CrownClosure == 0.0)
                return 0.0; //no evaporation through the snow cover...except for forest

            double actualEvaporation = potEvaporation;
            
            if(soilMoisture < _LP*_FC) //evaporation is supply limited
                actualEvaporation = potEvaporation * (soilMoisture / (_LP * _FC));
            
            return Math.Min(soilMoisture, actualEvaporation);
        }

#endregion


    }
}
